function ana_out = lmj_modeanalysis(funcmod, pertmat, data, in, varargin)
% XIPLOT_ALL
% Given multiple group of parameter values
% evaluate the model and plot autocovriance, variance decomposition, IRF,
% impulse responses.
% apply Kalman Filter to data and plot states, errors and
% state-decomposition
%
% INPUT:
%       "FUNCMOD":      function handle of the model
%       "PERTMAT":      [NPar NPert] matrix of parameters.
%       "DATA":         [NOBS NSER] matrix of data
%       "IN":           Structure that contains other necessary input
%          .STATE_SEL:  Cell array of names of selected state
%          .STANAMES:   names of all STATES
%          .SHONAMES:   names of all shocks
%          .SAVEPATH:   name of the save path
%          .PARNAMES:   name of parameters
%          .DATEVEC:    vector of dates
%          .MODELNAMES: names of the models
%          .FLAG_BLOCK: Whether to put
% OUTPUT:  LOTS and LOTS of graphs
%
%
% Dependencies (mfiles and matlab functions)
%   (major functions)
%       kfilter_all.m
%       momentszselstate.m
%       plot_multiple3d.m
%       accov.m
%       priorposterior_acout.m
%       plot_decompositions.m
%       multimod_acout.m (subfunction)
%
%   (tools) strdate, cr_dir, cellposition, dispaj, suptitle, ps2pdf, vline
%
% For: Alejandro Justiniano
% Research Assistant: Xi Luo
% Created: 02/22/2010

cucd = pwd;
% timestamp = strdate;

% Parse INPUTS
state_sel = in.state_sel;
stanames =  in.stanames;
shonames =  in.shonames;
mkdir(in.savepath, in.savefolder);
pert_savepath = [in.savepath, '\',in.savefolder];
cd(pert_savepath); mkdir('SinglePage'); cd(cucd);

[e, in] = ch_field(in, 'flag_plotst', 1);
[e, in] = ch_field(in, 'flag_plotss', 0);
[e, in] = ch_field(in, 'flag_irfbysho', 1);

flag_plotst = in.flag_plotst; 
flag_plotss = in.flag_plotss; 
flag_irfbysho = in.flag_irfbysho; 

parnames =  in.parnames;
datevec = in.dates;
modelnames = in.modelnames;
flag_block = in.flag_block;
[numpar, numpert] = size(pertmat);
state_pos = cellposition(state_sel, stanames);

% Solve for each PARAMETER.
numst = length(stanames);
numshock = length(shonames);
pertv_ok = zeros(numpert, 1);
GGmat = zeros(numst, numst, numpert);
CONSmat = zeros(numst, 1, numpert);
RRmat = zeros(numst, numshock, numpert);
SDXmat = zeros(numshock, numshock, numpert);
eumat = zeros(2, 1, numpert);
ZZmat = zeros(size(data,2),numst,numpert ); 

for ii = 1:numpert
    [GGtemp,RRtemp,CONStemp,eutemp,SDXtemp,ZZtemp]=feval(funcmod,pertmat(:, ii), varargin{:});
    if isequal(eutemp,[1;1])
        dispaj('pertvec Number ',ii,'  will be plotted');
        pertv_ok(ii)=1;
        GGmat(:, :, ii) = GGtemp;
        RRmat(:, :, ii) = RRtemp;
        CONSmat(:, :, ii) = CONStemp;
        eumat(:, :, ii) = eutemp;
        SDXmat(:, :, ii) = SDXtemp;
        ZZmat(:, :, ii) = ZZtemp;
    else
        dispaj('pertvec Number ',ii,'  cannot be plotted');
        pertv_ok(ii)=0;
    end
end
clear GGtemp RRtemp CONStemp eutemp SDXtemp ZZtemp ii

posok = find(pertv_ok == 1);
if length(posok(:)) < 1
    error('No parameter in PERTMAT delivers unique and stable solutions')
end

% Taking care of parameters that does not yield
% the stable and unique solutions.
pertmat_ok = pertmat(:, posok); % matrices of parameters that yield solution
GGmat_ok = GGmat(:, :, posok);
RRmat_ok = RRmat(:, :, posok);
CONSmat_ok = CONSmat(:, :, posok);
eumat_ok = eumat(:, :, posok);
SDXmat_ok = SDXmat(:, :, posok);
ZZmat_ok = ZZmat(:, :, posok);
numpert_ok = size(pertmat_ok, 2);
modelnames_ok = modelnames(posok);
dispaj(numpert_ok, 'out of', numpert, 'starting parameters deliver steady state');

% Variables such as "blabla_ok" are relevant in our analysis because when
% we perturb one element of a single parameter, we are not always sure if
% the newly constructed parameter will deliver stable solution

% when we load xmode from different models, we are more confident that
% parameters in PERTMAT will deliver unique and stable solutions.


%% Apply Kalman Filther/Smoother to parameters
clear stt_mat etamat_mat countmat_mat
[nobs] = size(data, 1);
stt_mat = zeros(nobs, numst, numpert_ok);
etamat_mat = zeros(nobs, numshock, numpert_ok);
countmat_mat = zeros(nobs, numst, numshock, numpert_ok);
for ii = 1:numpert_ok
    [stt,etamat,smooth_st,yfor,vstar,logLnc,countmat,countobs,countszero]= ...
        lmj_kfilterall(funcmod,pertmat_ok(:, ii),data, 1, [], varargin{:}); % flag_mode,ins);
    stt_mat(:, :, ii) = stt;
    etamat_mat(:, :, ii) = etamat;
    countmat_mat(:, :, :, ii) = countmat;
end


%% Starting the BIG LOOP HERE
ZZtemp = ZZmat(:, :, 1);
[posz, junk] = find(ZZtemp' ~= 0);
clear ZZtemp junk

sernames = stanames;
if ~exist('znames', 'var')
    znames = sernames(posz);
end

%%


clear stdthz_mat stdths_mat acthz_mat acths_mat covthz_mat covths_mat
clear vdecstaz_mat vdecstas_mat irfz_mat irfs_mat vdechz_mat vdechs_mat
clear stdsimz_mat stdsims_mat acsimz_mat acsims_mat covsimz_mat covsims_mat

for ii = 1:numpert_ok
    dispaj('begin parameter', ii);
    GG = GGmat_ok(:, :, ii);
    RR = RRmat_ok(:, :, ii);
    CONS = CONSmat_ok(:, :, ii);
    SDX = SDXmat_ok(:, :, ii);
    ZZ = ZZmat_ok(:, :, ii);
    out = momentszselstate(GG,CONS,RR,SDX,ZZ,in,state_pos);
    
    % STORing in MATRICES
    stdthz_mat(:, :, ii) = out.stdthz;
    stdths_mat(:, :, ii) = out.stdths;
    acthz_mat(:, :, :, ii) = out.acthz;
    acths_mat(:, :, :, ii) = out.acths;
    covthz_mat(:, :, :, ii) = out.covthz;
    covths_mat(:, :, :, ii) = out.covths;
    vdecstaz_mat(:, :, ii) = out.vdecstaz;
    vdecstas_mat(:, :, ii) = out.vdecstas;
    irfz_mat(:, :, :, ii) = out.irfz;
    irfs_mat(:, :, :, ii) = out.irfs;
    vdechz_mat(:, :, :, ii) = out.vdechz;
    vdechs_mat(:, :, :, ii) = out.vdechs;
    stdsimz_mat(:, :, ii) = out.stdsimz;
    stdsims_mat(:, :, ii) = out.stdsims;
    acsimz_mat(:, :, :, :, ii) = out.acsimz;
    acsims_mat(:, :, :, :, ii) = out.acsims;
    covsimz_mat(:, :, :, :, ii) = out.covsimz;
    covsims_mat(:, :, :, :, ii) = out.covsims;
end
clear GG RR CONS SDX ZZ ii

ana_out.stdsimz_mat = stdsimz_mat; 

nz = size(ZZmat_ok, 1);
ns = length(state_pos);


%% ========================================================================
%  START Plotting stuff!
% =========================================================================
%% Compare Volatilities
disp(' ')
disp('______________________________________________')
disp('Compare Volatilities of Observables with Real Data');
gvec=zeros(nz,1);
gnames=cell(nz,1);
ghand = figure;
for ii=1:nz;
    for jj=1:numpert_ok
        temp_base=std(data(:,ii));
        if numpert_ok ~=1
            subplot(nz,numpert_ok,jj + (ii-1)*numpert_ok);
            
        else
            subplot(3,ceil(nz/3),jj + (ii-1)*numpert_ok);
            title(znames{ii});
        end
        temp = stdsimz_mat(ii, :, jj);
        hist(temp,20);
        if numpert_ok ~= 1
            if ii == 1; title(modelnames_ok{jj}); end
            ylabel(znames{ii});
        else
            title(znames{ii})
        end
        hand=vline(temp_base, 'r-');
        set(hand,'LineWidth',4);
        
        
    end
    %
    %     gvec(ii)=gcf;
    %     gnames(ii)={['Z STD ',znames{ii}]};
end

if numpert_ok ~= 1;
    suptitle('Z STD ',14);
else
    suptitle(['Z STD ', modelnames_ok{1}], 14);
end

subtitulo(pert_savepath);

% set(ghand,'PaperOrientation','Portrait');
set(ghand,'PaperPosition',[0.05 0.03 8.5 11]);

gname = 'Z STD';
cd(pert_savepath);
% print pages one by one
print(ghand, '-dpdf', gname);

cd(cucd);
close all;


%% IRFS
disp(' ')
disp('______________________________________________')
disp('Plotting IRF of Observables')
    nc = 3; 
    ss.subdesign(2)=nc; 
    znames_full = cell(nz, 1); 
    for ii = 1:nz; 
        if in.flag_addz(ii) == 0; 
            znames_full{ii} = znames{ii};
        elseif in.flag_addz(ii) == 1; 
            znames_full{ii} = ['lev ', znames{ii}];
        end 
    end
        
    
    
    
if flag_irfbysho == 1;
    ss.subdesign(1)=ceil(nz/nc);
    ss.rnames =znames_full;  
    ss.cnames = {''};
    ss.bnames = modelnames_ok;
    gvec = zeros(length(shonames), 1);
    gnames = cell(length(shonames), 1); 
    for ii = 1:length(shonames);
        ss.suptitle = ['Z IRF: ', shonames{ii}]; 
        gvec(ii) = plot_multiple3d(irfz_mat(:,ii,:,:),ss); 
        subtitulo(pert_savepath);
        gnames(ii)={['Z IRF ',shonames{ii}]};
    end 
   
    
else
    ss.subdesign(1)=ceil(length(shonames)/nc);
    ss.rnames=znames;
    ss.cnames=shonames;
    ss.bnames=modelnames_ok;
    ss.legfont=11;
    ss.suptitfont=14;
    ss.titfont=12;
    
    gvec=zeros(nz,1);
    gnames=cell(nz,1);
    for ii=1:nz;
        if in.flag_addz(ii) == 0
            ss.suptitle=['Z Impulse response   ',znames{ii}];
        elseif in.flag_addz(ii) == 1
            ss.suptitle=['Z Impulse response: level of ',znames{ii}];
        end
        % ss.suptitle=['Z Impulse response ',znames{ii}];
        ss.rnames={''};
        gvec(ii)=plot_multiple3d(irfz_mat(ii,:,:,:),ss);
        subtitulo(pert_savepath);
        
        gnames(ii)={['Z IRF ',znames{ii}]};
    end
    
    
end

cd(pert_savepath)
    cd SinglePage
    for ii=1:length(gvec);
        print(gvec(ii),'-dpdf',gnames{ii});
    end
    cd(pert_savepath)
    
    sumname = 'Sum Z IRF';
    print(gvec(1), '-dpsc', sumname);
    for ii = 2:nz;
        print(gvec(ii), '-dpsc', '-append', sumname);
    end
    ps2pdf('psfile', [sumname,'.ps'], 'pdffile', [sumname, '.pdf'])
    
    cd(cucd);
    close all;

disp(' ')
disp('______________________________________________')
disp('Plotting IRF of Selected States')
ss.subdesign(1)=5;
ss.subdesign(2)=3;
ss.rnames=sernames();
ss.cnames=shonames;
ss.bnames=modelnames_ok;
ss.legfont=11;
ss.suptitfont=14;
ss.titfont=12;

gvec=zeros(ns,1);
gnames=cell(ns,1);
for ii=1:ns;
    if in.flag_adds(ii) == 0
        ss.suptitle=['S Impulse response   ',sernames{state_pos(ii)}];
    elseif in.flag_adds(ii) == 1
        ss.suptitle=['S Impulse response: level of ',sernames{state_pos(ii)}];
    end
    ss.rnames={''};
    gvec(ii)=plot_multiple3d(irfs_mat(ii,:,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(ii)={['S IRF ',sernames{state_pos(ii)}]};
end

cd(pert_savepath)
cd SinglePage
for ii=1:length(gvec);
    print(gvec(ii),'-dpdf',gnames{ii});
end
cd(pert_savepath)

sumname = 'Sum S IRF';
print(gvec(1), '-dpsc', sumname);
for ii = 2:ns;
    print(gvec(ii), '-dpsc', '-append', sumname);
end
ps2pdf('psfile', [sumname,'.ps'], 'pdffile', [sumname, '.pdf'])

cd(cucd);
close all;

%%  VDEC
disp(' ')
disp('_______________________________________________')
disp('Variance Decomposition for Observables')

gvec=zeros(nz,1);
gnames=cell(nz,1);

for ii=1:nz;
    ss.suptitle=['Z Variance decomposition ' znames{ii}];
    ss.rnames={''};
    gvec(ii)=plot_multiple3d(vdechz_mat(ii,:,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(ii)={['Z VDEC ',znames{ii}]};
end

cd(pert_savepath)
cd SinglePage
for ii=1:length(gvec);
    print(gvec(ii),'-dpdf',gnames{ii});
end
cd(pert_savepath)

sumname = 'Sum Z VDEC';
print(gvec(1), '-dpsc', sumname);
for ii = 2:nz;
    print(gvec(ii), '-dpsc', '-append', sumname);
end
ps2pdf('psfile', [sumname,'.ps'], 'pdffile', [sumname, '.pdf'])
cd(cucd);
close all;


%% Variance Decomposition for Selected States

disp(' ')
disp('_______________________________________________')
disp('Variance Decomposition for Selected States')

gvec=zeros(ns,1);
gnames=cell(ns,1);

for ii=1:ns;
    ss.suptitle=['S Variance decomposition ' sernames{state_pos(ii)}];
    ss.rnames={''};
    gvec(ii)=plot_multiple3d(vdechs_mat(ii,:,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(ii)={['S VDEC ',sernames{state_pos(ii)}]};
end
cd(pert_savepath)
cd SinglePage
for ii=1:length(gvec);
    print(gvec(ii),'-dpdf',gnames{ii});
end
cd(pert_savepath)

sumname = 'Sum S VDEC';
print(gvec(1), '-dpsc', sumname);
for ii = 2:ns;
    print(gvec(ii), '-dpsc', '-append', sumname);
end
ps2pdf('psfile', [sumname,'.ps'], 'pdffile', [sumname, '.pdf'])
cd(cucd);
close all;


%% Autocorrelation

disp(' ')
disp('________________________________________________')
disp('Autocorrelation for Observables')

ac_data=accov(data,(0:in.nper));
clear med_acsimz_mat;
for ii = 1:numpert_ok+1;
    if ii == 1;
        med_acsimz_mat(:,:,:,ii)=ac_data;
    else
        med_acsimz_mat(:,:,:,ii)=median(acsimz_mat(:, :, :, :, ii-1),4);
    end
end
%med_acsimz_mat(:,:,:,numpert_ok + 1)=ac_data;
% Partiotion into four blocks

if flag_block
    gvec=zeros(4,1);
    gnames=cell(4,1);
    block1=(1:5);
    block2=(6:nz);
else
    gvec = zeros(1, 1);
    gnames = cell(1, 1);
    block1 = (1:nz);
end

try
    ss=rmfield(ss,'subdesign');
catch
    disp('subdesign already removed')
end

ss.rnames=sernames(posz(block1));
ss.bnames=[{'data'}; modelnames_ok(:)];
ss.legfont=9;
ss.titfont= 11;
if flag_block
    ss.suptitle='Obs AC Median Block 1 1';
else
    ss.suptitle='Obs AC Median Block all';
end

ss.cnames=sernames(posz(block1));
gvec(1)=plot_multiple3d(med_acsimz_mat(block1,block1,:,:),ss);
subtitulo(pert_savepath);

if flag_block
    gnames(1)={'Obs AC MED 1 1'};
else
    gnames(1)={'Obs AC MED all'};
end

if flag_block
    ss.suptitle='Obs AC Median Block 1 2';
    ss.cnames=sernames(posz(block2));
    gvec(2)=plot_multiple3d(med_acsimz_mat(block1,block2,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(2)={'Obs AC MED 1 2'};
    
    ss.suptitle='Obs AC Median Block 2 2';
    ss.rnames=sernames(posz(block2));
    ss.cnames=sernames(posz(block2));
    gvec(3)=plot_multiple3d(med_acsimz_mat(block2,block2,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(3)={'Obs AC MED 2 2'};
    
    ss.suptitle='Obs AC Median Block 2 1';
    ss.rnames=sernames(posz(block2));
    ss.cnames=sernames(posz(block1));
    gvec(4)=plot_multiple3d(med_acsimz_mat(block2,block1,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(4)={'Obs AC MED 2 1'};
end

cd(pert_savepath)
if ~flag_block
    print(gvec(1),'-dpdf',gnames{ii});
else
    gname = 'Obs AC Median by Block';
    print(gvec(1), '-dpsc', gname);
    for ii = 2:length(gvec);
        print(gvec(ii), '-dpsc', '-append', gname);
    end
    ps2pdf('psfile', [gname,'.ps'], 'pdffile', [gname, '.pdf'])
    
end
cd(cucd);
close all;


%% Autocorrelation for Selected States

disp(' ')
disp('________________________________________________')
disp('Autocorrelation for Selected States')

ac_data=accov(data,(0:in.nper));
clear med_acsims_mat;
for ii = 1:numpert_ok;
    med_acsims_mat(:,:,:,ii)=median(acsims_mat(:, :, :, :, ii),4);
end

% Partiotion into four blocks... or not

if flag_block
    gvec=zeros(4,1);
    gnames=cell(4,1);
    block1=(1:3);
    block2=(4:ns);
else
    gvec = zeros(1, 1);
    gnames = cell(1, 1);
    block1 = (1:ns);
end

try
    ss=rmfield(ss,'subdesign');
catch
    disp('subdesign already removed')
end

ss.rnames=sernames(state_pos(block1));
ss.bnames=modelnames_ok(:);
ss.legfont=9;
ss.titfont= 11;
if flag_block
    ss.suptitle='Sel States AC Median Block 1 1';
else
    ss.suptitle='Sel States AC Median Block all';
end
ss.cnames=sernames(state_pos(block1));
gvec(1)=plot_multiple3d(med_acsims_mat(block1,block1,:,:),ss);
subtitulo(pert_savepath);

if flag_block
    gnames(1)={'States AC MED 1 1'};
else
    gnames(1)={'States AC MED all'};
end

if flag_block
    ss.suptitle='Sel States AC Median Block 1 2';
    ss.cnames=sernames(state_pos(block2));
    gvec(2)=plot_multiple3d(med_acsims_mat(block1,block2,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(2)={'States AC MED 1 2'};
    
    ss.suptitle='Sel States AC Median Block 2 2';
    ss.rnames=sernames(state_pos(block2));
    ss.cnames=sernames(state_pos(block2));
    gvec(3)=plot_multiple3d(med_acsims_mat(block2,block2,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(3)={'States AC MED 2 2'};
    
    ss.suptitle='Sel States AC Median Block 2 1';
    ss.rnames=sernames(state_pos(block2));
    ss.cnames=sernames(state_pos(block1));
    gvec(4)=plot_multiple3d(med_acsims_mat(block2,block1,:,:),ss);
    subtitulo(pert_savepath);
    
    gnames(4)={'States AC MED 2 1'};
end

cd(pert_savepath)
if ~flag_block
    print(gvec(1),'-dpdf',gnames{ii});
else
    gname = 'States AC Median by Block';
    print(gvec(1), '-dpsc', gname);
    for ii = 2:length(gvec);
        print(gvec(ii), '-dpsc', '-append', gname);
    end
    ps2pdf('psfile', [gname,'.ps'], 'pdffile', [gname, '.pdf'])
    
end
cd(cucd);
close all;

%% Autocorrelation Full for each

gvec=zeros(numpert_ok + 1,1);
gnames=cell(numpert_ok+1,1);
for ii = 1:numpert_ok
    gvec(ii) = priorposterior_acout(acsimz_mat(:, :, :, :, ii),ac_data, sernames(posz), ...
        [], (1:nz), (1:nz), [0.05, 0.95]);
    
    suptitle(modelnames_ok{ii});
    subtitulo(pert_savepath);
    
    gnames{ii} = ['par', num2str(ii), ' AC all Z'];
    hand = gcf;
    set(hand,'PaperOrientation','Landscape');
    set(hand,'PaperPosition',[0.2 0.1 10.5 8.3]);
end

cd(pert_savepath)
cd SinglePage
for ii = 1:numpert_ok
    print(gvec(ii),'-dpdf',gnames{ii});
end
cd(cucd)


gvec(numpert_ok + 1) = multimod_acout(acsimz_mat,ac_data, ...
    sernames(posz), [modelnames_ok, 'data'], (1:nz), (1:nz));

suptitle('comparison of AC accross modes');
subtitulo(pert_savepath);

gnames{numpert_ok + 1} = 'compare modes AC all Z';
hand = gcf;
set(hand,'PaperOrientation','Landscape');
set(hand,'PaperPosition',[0.2 0.1 10.5 8.3]);
cd(pert_savepath)
cd SinglePage
print(gvec(numpert_ok + 1),'-dpdf',gnames{numpert_ok + 1});
cd(pert_savepath)

sumname = 'Sum AC all Z';
print(gvec(1), '-dpsc', sumname);
for ii = 2:numpert_ok+1;
    print(gvec(ii), '-dpsc', '-append', sumname);
end

ps2pdf('psfile', [sumname,'.ps'], 'pdffile', [sumname, '.pdf'])

cd(cucd)
close all;

%% STATE for different modes
if flag_plotst == 1
    stt_mat(:, :, ii) = stt;
    disp('________________________________________________')
    disp('Selected States for different modes')
    sel_sttmat = stt_mat(:, state_pos, :);
    num_selst = length(state_pos);
    hand = figure;
    nc = 3; 
    nr = ceil(num_selst/nc); 
    colorvec = {'b', 'g', 'r', 'c', 'b:', 'g:', 'r:', 'c:'};
    for ii = 1:num_selst
        subplot(nr, nc, ii)
        for jj = 1:numpert_ok
            plot(sel_sttmat(:, ii, jj), 'color', colorvec{jj});
            if jj == 1; hold on; end
            if jj == numpert_ok; hold off; end
        end
        title(state_sel{ii});
    end
    subplot(nr, nc, nc*nr)
    for ii = 1:numpert_ok
        plot(ones(3, 1), 'color', colorvec{ii});
        if ii == 1; hold on; end
    end
    legend(modelnames_ok)
    suptitle('State for Different Models')
    subtitulo(pert_savepath);
    
    set(hand,'PaperOrientation','Landscape');
    set(hand,'PaperPosition',[0.2 0.1 10.5 8.3]);
    gname = 'KF State';
    cd(pert_savepath); print(hand,'-dpdf',gname); cd(cucd);
    
    %% Decomposition of COUNTMAT
    gvec = zeros(numpert_ok*num_selst, 1);
    sel_countmatmat = countmat_mat(:, state_pos, :, :);
    gcounter = 0;
    for ii = 1:numpert_ok;
        for jj = 1:num_selst
            tit = ['STATE: ', state_sel{jj}, ';  MODEL: ', modelnames_ok{ii}];
            gcounter = gcounter + 1;
            gvec(gcounter) =plot_decompositions(squeeze(sel_sttmat(:, jj, ii)), ...
                squeeze(sel_countmatmat(:, jj, :, ii)) ,state_sel{jj},shonames,datevec,tit,0,0);
            gname = ['STDECOMP_',   state_sel{jj}, '_model', num2str(ii)];
            subtitulo(pert_savepath);
            
            cd(pert_savepath); cd SinglePage; print(gvec(gcounter),'-dpdf',gname); cd(cucd);
        end
    end
    
    cd(pert_savepath)
    sumname = 'ST_DecompAll';
    print(gvec(1), '-dpsc', sumname);
    for ii = 2:numpert_ok*num_selst;
        print(gvec(ii), '-dpsc', '-append', sumname);
    end
    ps2pdf('psfile', [sumname,'.ps'], 'pdffile', [sumname, '.pdf'])
    cd(cucd)
    close all
    
    %% Plot ERROR terms
    etamat_mat(:, :, ii) = etamat;
    disp('__________________________________________')
    disp('Error Terms')
    num_shocks = length(shonames);
    hand = figure;
    nc = 3; 
    nr = ceil(num_shocks/nc); 
    for ii = 1:num_shocks
        subplot(nr, nc, ii)
        for jj = 1:numpert
            plot(etamat_mat(:, ii, jj), 'color', colorvec{jj})
            if jj == 1; hold on; end
            if jj == numpert; hold off; end
        end
        title(shonames{ii})
        if ii == num_shocks; legend(modelnames_ok, 'Location','Best'); end
    end
    suptitle('Error Terms accross models ')
    subtitulo(pert_savepath);
    
    gname = 'KF errors';
    set(hand,'PaperOrientation','Landscape');
    set(hand,'PaperPosition',[0.2 0.1 10.5 8.3]);
    cd(pert_savepath); print(hand,'-dpdf',gname); cd(cucd);
end

%% Plot Steady state
 


if flag_plotss == 1 && numpert_ok >= 2;
    disp('________________________________________________')
    disp('Selected Steady States for different modes')
    ssreport = in.ssreport; 
    ssnames  = in.ssnames; 
    pospert =  in.pospert; 
    pertvec_ok = pertmat_ok(pospert, :);
    pertvec_temp = sort([pertmat(pospert, 1), min(pertvec_ok),max(pertvec_ok)]); 
    pertbnd = [pertvec_temp(1), pertvec_temp(end)]; 
    pertvec_long = sort(unique([pertmat(pospert, 1), linspace(pertbnd(1), pertbnd(2), 10)])); 
    numpert_long = length(pertvec_long); 
    pertmat_long = repmat(pertmat(:, 1), 1, numpert_long); 
    pertmat_long(pospert, :) = pertvec_long; 
    ssvec_mat = nan(length(ssnames), numpert_long); 
    for ii = 1:numpert_long
        [GGtemp,RRtemp,CONStemp,eutemp,SDXtemp,ZZtemp, initss, ssvec]=feval(funcmod,pertmat_long(:, ii), varargin{:});
        clear GGtemp RRtemp CONStemp SDXtemp ZZtemp initss
        if isequal(eutemp,[1;1])
            dispaj('pertvec Number ',ii,'  will be plotted');
            ssvec_mat(:, ii) = ssvec;
        else
            dispaj('pertvec Number ',ii,'  cannot be plotted');
        end
        clear eutemp
    end
    
    num_ssreport = length(ssreport); 
    pos_ssreport = cellposition(ssreport, ssnames); 
    ssvec_mat_report = ssvec_mat(pos_ssreport, :);
   
    xtrue = pertmat(pospert, 1); 
    vlinestr = modelnames{1}; 
    nc = 2; 
    nr = ceil(num_ssreport/2); 
    ghand = figure; 
    
    for ii = 1:num_ssreport; 
        subplot(nr, nc, ii)
        y = ssvec_mat_report(ii, :); 
        plot(pertvec_long, y, 'LineWidth',4)
        vline(xtrue, 'r-', vlinestr);
        title(ssreport{ii})
        xlabel(parnames{pospert})
    end 
    suptitle(['Sel Steady States: ', parnames{pospert}])
    subtitulo(pert_savepath);
    set(ghand,'PaperPosition',[0.05 0.03 8.5 11]);
    gname = 'Sel Steady States';
    cd(pert_savepath);
    print(ghand, '-dpdf', gname);
    cd(cucd);     
    
end


%% 


close all
cd(pert_savepath)
delete *.ps
cd(cucd)
winopen(pert_savepath)


%% =====================================================================
% Begin NESTED FUNCTIONS
%  =====================================================================
    function hand=multimod_acout(matdr,matdata,repnames,legcell,blockr,blockc)
        % =========================================================================
        % This is like wedges_outcorr but for the project with Bruce
        % Key difference is ability to determine bands and also which blocks to
        % plot against one another
        
        % priorposterior_acout(matdr,matdata,repnames,blockr,blockc,pvec) % wedges_outcorr
        % MATDR         [NR NC NP NDR NPERT]
        % MATDATA       [NR NC NP    ]
        % Code assumed NR=NC, can esily be generalized
        % REPNAMES      [NR 1] Cell of names
        % BLOCKR        Index of NR that will make the rows
        % BLOCKC        Index of NC that will make the columns
        % Xi Luo 02/19/2010,
        % revised from AJ's function priorposterior_acout.m (October 17 2006)
        %
        % =========================================================================
        dimvec=size(matdr);
        if length(dimvec) < 4;
            error('MATDR must be 4D');
        end
        [nr nc ncov ndr npert]=size(matdr);
        xaxis=(0:ncov-1)';
        
        if ~isequal( size(matdata), [nr nc ncov] );
            error('MATDATA and MATDR do not match')
        end
        
        nrow=length(blockr);
        ncol=length(blockc);
        
        titfont=9;
        set(0,'DefaultAxesFontSize',8);
        % Transform the covariance matrix into a correlation matrix
        
        matdr=sort(matdr,4);
        matmed=median(matdr,4);
        %matmean=mean(matdr,4);
        
        figure;
        gind=0;
        for ii_sub=1:nrow
            for jj_sub=1:ncol
                gind=gind+1;
                subplot(nrow,ncol,gind);
                for kk = 1:npert
                    hand=plot(xaxis,squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, kk)),':*','LineWidth',2.5);hold on;
                    set(hand,'Color',(1 - 0.9*kk/npert)*[0.7, 0.3, 0.6],'MarkerSize',4)
                end
                hand=plot(xaxis,squeeze(matdata(blockr(ii_sub),blockc(jj_sub),:)),'LineWidth',2.5);
                set(hand,'Color',[0.1 0.1 0.1]);
                hold off;
                title([repnames{blockr(ii_sub)},'_t , ',repnames{blockc(jj_sub)},'_t_-_k'],'FontSize',titfont);
                box off;
                if npert > 1
                    ymax=max(max([squeeze(matdata(blockr(ii_sub),blockc(jj_sub),:))',squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, 1))', ...
                        squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, 2))']));
                    ymin=min(min([squeeze(matdata(blockr(ii_sub),blockc(jj_sub),:))',squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, 1))', ...
                        squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, 2))']));
                elseif npert == 1;
                    ymax=max(max([squeeze(matdata(blockr(ii_sub),blockc(jj_sub),:))',squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, 1))']));
                    ymin=min(min([squeeze(matdata(blockr(ii_sub),blockc(jj_sub),:))',squeeze(matmed(blockr(ii_sub),blockc(jj_sub),:, 1))']));
                end
                
                ylim([ymin-0.02 ymax+0.02]);
                
                
                if ncov < 8
                    set(gca,'XTick',(0:ncov-1));
                else
                    set(gca,'XTick',(0:round(ncov/2):ncov-1));
                end
                
                if ii < nrow
                    set(gca,'XTickLabel','')
                end
                
                xlim([0 (ncov-1)+0.1]);
                if ii==nrow && jj== ncol  && ~isempty(legcell);
                    hand=legend(legcell,'Location','SouthOutside');
                    set(hand,'box','off');
                    set(hand,'FontSize',7,'FontWeight','Bold');
                end
            end
        end
        
        hand=gcf;
        set(hand,'PaperOrientation','Landscape');
        set(hand,'PaperPosition',[0.1 0.1 10.92 7.8])
    end

end